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^-^ . A discrete-time model of reacting evolving fields, transported by a bidi- 



mensional chaotic fluid flow, is studied. Our approach is based on the use of a 
Lagrangian scheme where fluid particles are advected by a 2d symplectic map 



K> , possibly yielding Lagrangian chaos. Each fluid particle carries concentrations 

of active substances which evolve according to its own reaction dynamics. 
This evolution is also modeled in terms of maps. Motivated by the question, 
of relevance in marine ecology, of how a localized distribution of nutrients or 
preys affects the spatial structure of predators transported by a fluid flow, we 
study a specific model in which the population dynamics is given by a logistic 
map with space-dependent coefficient, and advection is given by the standard 
map. Fractal and random patterns in the Eulerian spatial concentration of 
predators are obtained under different conditions. Exploiting the analogies of 



this coupled-map (advection plus reaction) system with a random map, some 
features of these patterns are discussed. 



The spatial structure of passive fields transported by a fluid flow is an im- 
portant question in fluid and nonlinear dynamics. Many important situations 
require in addition to take into account chemical or biological interactions be- 
tween the substances transported by the flow. In particular, our motivation 
comes from the general question on the mechanisms for plankton inhomogeneity 
in the ocean, and more specially for the spatial patterns that would be reached 
by plankton or marine hervibores grazing from a localized source of nutrients. 
We present a general methodology, based in coupled discrete-time dynamical 
systems describing reaction or population dynamics and advection in the La- 
grangian framework, and apply it to a model of logistic population dynamics in 
the presence of a localized source of nutrients and chaotic advection. Fractal 
and random spatial features develop in the concentration patterns, which can 
be understood in terms of an analogy with random maps. 

I. INTRODUCTION 

Patchiness, or the uneven distribution of substances of organisms, is ubiquitously ob- 
served in the ocean [0-§]. In complex situations such as the marine ecosystems, characterized 
by the interplay of population dynamics and an ambient fluid motion which may differently 
affect individual populations, the question of how a localized availability of nutrients (or 
preys, or, in chemical terms, activators) may affect the distribution of primary producers 
(or predators, or inhibitors, respectively; in the following we refer to preys and predators) 
is a crucial and challenging problem. 

According to 0, the possible causes originating patchiness in marine ecosystems may 
be grouped in different categories: fluid motion, biological growth coupled with dispersive 
processes, less ubiquitous mechanisms like swarming, vertical migration, and others. In 
addition to patchiness, localized availability of food may occur in correspondence of localized 
sub-ecosystems, such as Posidonia Oceanica beds (for a review see 0): they display a quite 



complex structure in which the most important features arise from direct consumption of 
the plant and epiphyte-herbivore interactions, although a portion of the trophic chain is 
based on suspended matter |Q. 

A comprehensive description of such processes leads to the study of the so called 
advection-reaction-diffusion equations. These are partial differential equations of the type: 



dt 



(v(x,t) ■ V)Q = Ri{Cu ...,CAr,x,t) + AAQ. (1) 



Where Ci{'}c,t) {i = 1, ..., A^) is the concentration of the i-th reactive (in biological or chemical 
terms) species or substance, and the functions Ri describe the reaction, or the population, 
intrinsic dynamics. The possible explicit spatial and temporal dependence may model the 
influence of temporal and spatial inhomogeneities in food, temperature, etc. The term 
(v(x, t) • V)Cj represents advection by a given solenoidal (i.e. incompressible V ■ v = 0) 
velocity field v(x, t). Finally, the term DiACi describes diffusion of the i-th species or 
substance with diffusivity Di. In writing (|1]) we are assuming that the evolution of the 
advected concentrations does not affect that of the underlying flow v(x, t). This is definitely 
reasonable at the scales we are interested in, even though it is worth mentioning that at much 
smaller scales the presence of organisms may affect the rheological properties of seawater 

ii- 

In some situations it would be necessary to consider a different velocity field for each 
of the A^ concentrations Cj. This may happen, for instance, if different organisms live at 
different mean depths, leading to differences in the experienced flow. In such cases one has 
to replace Eq. (|I]) by 

g^' ' + (v.(x, t) ■ V)C, = i?,(Ci, ..., Cn, X, t) + AAC,. (2) 

In this paper, however, we will only consider the case given by Eq. (|I|) in which the same 
velocity field advects all the substances. 

Introducing the Lagrangian time derivative -^ = ^ + v ■ V, Eq. (^ can be written in the 
form 



^^^ = i?,(Ci,...,C^,t) + AAC,. (3) 

If diffusion is neglected, it is simple to write the Ci(x, t) in terms of the solutions of the 
Lagrangian evolution equation 

^ = v(x(t),t), (4) 

and the reactive evolution equation 

^^ = i?,(Ci,...,C^,x(t),t). (5) 

This set of coupled ordinary differential equations describe the advection-reaction process 
in a Lagrangian frame: fluid particles move according to (^), and reactions among the Cj's 
occur inside each fluid particle, as expressed by (|^), where Ci{t) are the concentrations at a 
particular fluid particle (the one at x{t) at time t, i.e. Ci(t) = Cj(x(t),t)). Denoting by S* 
and L* the formal solutions of @ and (H) respectively (i.e., x(t) = S'*x(0) and C(t) = L*C(0), 
with C = (Ci, C2, ..., Cn)) we can write the solution of Eq. (|l]) in the form (with Di = 0) 

C(x,t) = L*C(^-*x,0). (6) 

The case Di ^ needs a more elaborated treatment. In addition to the time dependence, 
L* will have also an explicit space dependence if the i?j's have it. 

Obviously, the detailed understanding of the above class of partial differential equations 
constitute a formidable task. However, at this stage of development, we are just interested 
in the search for generic behaviors expected when few typical characteristics of the flow and 
of the population dynamics are considered. Thus, if, for example, we concentrate on flows of 
geophysical nature, horizontal motion turns out to be much more intense than vertical one 
as soon as one considers scales larger than a few kilometers. This justifies restricting in the 
following to incompressible twodimensional flows. A turbulent bidimensional flow would be 
a way to model the irregular advection process to which suspended matter is subjected in 
real oceans. There are however simpler classes of flows which share some basic character- 
istics with turbulence, but are much more accessible to analysis: Lagrangian chaotic flows 



PJ^. These are smooth velocity fields, with some simple time dependence in the Eulerian 
description, but which lead to chaotic trajectories of fiuid elements, with the associated 
stretching and folding, in the Lagrangian description. In this restricted framework, it is well 
known that even periodic time-dependence in two-dimensional incompressible fiows leads 
generically to chaotic motion of fiuid particles. 

Rather than integrating the full equations describing the continuous in time dynamics, 
and since our interest lies mainly in a qualitative characterization of the population system, 
we will resort in this paper to a discrete in time mapping-approach in terms of discrete- 
time dynamical systems. This approach is numerically very efficient, and has proven to 



be extremely productive to study the impact of chaotic advection on mixing ||8|-[T0|. The 
main idea is to mimic the advection and reaction processes in terms of maps that capture 
the main features of each aspect. Thus, since we will be looking at processes taking place 
in 2d incompressible fiows, the advective part of our model is naturally described by a 
two dimensional symplectic map. It is well known that Lagrangian motion in such systems 
is typically chaotic and with a rather rich behavior. In addition, the transported fiuid parcel 
contains concentrations of active chemical substances or biological specimens subjected to 
a specific dynamics that will be also modeled in terms of a map. Also, it is worth noticing 
that even though this is not done in this work, diffusion can be easily reincorporated into 
the model by averaging, after each iteration of the maps, the concentration of the different 
species contained in the fluid elements over a region of size / ~ \fD^, where Dj is the 
corresponding diffusivity and r is a characteristic time scale of the system. 

The general ideas and formalism sketched above will be made more concrete in the 
following, and applied to tackle our main problem: the influence of inhoniogeneities of the 
distribution of preys on that of predators. We will make use of some known results for random 
maps, and compare our discrete-time approach with results obtained in a continuous-time 



description of the problem [p!5| , p!6| . In Section II we will discuss the approach to population 
dynamics in terms of maps; in Section III our particular model is presented and compared 
with results for a random logistic map; the analogy helps in the interpretation of the resulting 
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spatial structure of predators, which is described in Section IV, and discussed in Section V. 

II. A DISCRETE-TIME APPROACH 

Let us now present the general idea of our approach for the analysis of the population 
dynamics in terms of maps. 

It is easy to understand that for time-periodic velocity fields, i.e. v(x, t) = v(x, t + T) 
where T is the period, Eq (^ can be described by a discrete-time dynamical system. The 
position :K{t+T) is univocally determined by x(t). In addition (because of the periodic veloc- 
ity field) the map x(t) —>■ x(t + T) cannot depend on t. Since a periodic time dependence is 
enough to induce Lagrangian chaos, and because of the above mathematical simplifications, 
we particularize our study to time-periodic velocity fields, for which we can write 

x(t + l)=F(x(t)), (7) 

Now, time is measured in units of the period T. If v is incompressible, the map (^ is 
volume (area in 2d) preserving, i.e., det(||^) = 1. In 2d the map (0) is symplectic, i.e. 
the discrete-time version of a Hamiltonian system. Usually, it is not simple at all to obtain 
F(x(t)) for a given v(x, t). However, one can directly write models for F which contain the 
qualitative features of the flow one is trying to model. 

In addition, the transported fluid parcel contains species subjected to their own popu- 
lation dynamics. Denoting the solution of (^ after one period of the flow (L^) by G, the 
evolution rule for the interacting concentrations C = (Ci, C2, ...C„) is expressed in terms of 
a map: 

C(t + l) = G(C(t)) (8) 

As before, G will carry additional explicit time and space dependencies if (|^) is not 
autonomous in space or time. The discrete-time version of Eq. (^ is: 

C(F(x),t + l) = G(C(x,t)). (9) 

In the following we particularize this general approach to a particular model. 



III. A PARTICULAR MODEL AND ITS RELATIONSHIP WITH A RANDOM 

LOGISTIC MAP 

The main interest of our study is to consider the problem of how the spatial structure of 
the prey spatial distribution may affect the one of the predators. In this section we study a 
particular model and present its analogies with the random map. This will enable us to use 
the already known properties of this random map to get a further insight into the influence 
of the distribution of preys on the predator patterns. We will considered a single-species 
population dynamics, i.e. the predator evolves for fixed prey distribution, and under the 
influence of the flow, but the distribution of the prey is a non-dynamic variable, in the sense 
that it is not transported by the flow and is maintained at fixed values undisturbed by the 
predator action. This is the simplest setting in which the effects of a localized source of 
nutrients on an advected predator will show up. 

The model is the following: the positions of the fluid parcels are advected by a standard 



map 1 11], i.e. a 2d symplectic map defined in the square of side 27r by 



x{t + 1) = (x(t) + K sin y{t + 1)) mod 2n (10) 

y{t + 1) = {y{t) + x{t)) mod 27r (11) 

It is not integrable for K ^ 0. As K increases chaotic regions occupy larger areas, and 
the original KAM tori (regular non-chaotic orbits) are successively destroyed. For K large 
enough the KAM tori occupy a very small region and practically the whole phase space is a 
unique chaotic region. 

The model is completed by stating the evolution rules for the predator-prey interactions. 
We denote with n(x) the stationary spatially inhomogenous distribution of preys, and with 
C(x, t) the concentration of predators in point x at time t. We take it to evolve in each 
fluid parcel according to the well known logistic map : C{t + 1) = G{C) = rC{t){l — C(t)), 
but with a growth rate parameter r determined by the presence of preys, i.e., r = fin. The 
complete evolution equation (y) is 
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C(F(x), t + 1) = /in(x)C(x, t)(l - C(x, t)). (12) 

The standard map has been written in the compact form (J^) with x = {x,y). 
We now introduce the particular form of the locahzed prey distribution n(x): 

\ ri iixe[n{l-p),n{l+p)], yy 
r = fin['x) = < (13) 

ro otherwise. 

with < p < 1. We are suggesting a striped spatial distribution of the preys (with strip 
width 2ttp), which basically represents (due to the 27r-periodicity of the flow) the simplest, 
space-periodic fashion to mimic a patchy distribution. A fraction p of system area has the 
value r = ri, and fraction 1 — p, the r = tq. The heuristic idea that will guide our analysis 
is that, if mixing provided by the advection map is strong enough, fluid parcels will visit 
regions with the different values of r in a stochastic way, so that the Lagrangian evolution of 
the concentrations will be well described by a random logistic map, i.e. a map of the form 

C{t+l)=atC{t)il-Cit)), (14) 

where the random variable at can take only two values 

{ro with probability 1 — p, 
(15) 
ri with probability p, 

and Qf+i is independent of any previous a^. 

The random map (p!^ has been studied in |T2| for ro = 1/2 and ri = 4. This corresponds 



to the situation in which for a value of r = ro the population dynamics is attracted by a 
fixed point, whereas chaotic population dynamics occurs for r = ri. The alternancy in time 
of these two tendencies gives rise to nontrivial behavior. Numerical (and some analytical) 



computations |jT2[ give the following results for ([l^ with ro = 1/2 and ri = 4, for different 
values of p: 



i) li p < Pi = 1/3 the Lyapunov exponent A = limTv^oo ;^ Z]j=o In joi (1 — 2C(t))| is 
negative, i.e. there is exponential convergence of two initially close sequences C{t) 



generated with the same sequence at but shghtly different initial conditions C(0). In 
this case, the sequences are attracted by C = 0. 

ii) If pi < p < p2 — 0.5 the Lyapunov exponent is negative again, but now the sequences 
do not converge to any fixed point. They wander in an irregular and seemingly chaotic 
manner (of the on-off interniittency type). The meaning of the negative value of 
A is that close initial values of C evolve, under the same sequence at, towards the 
same irregular trajectory. This is a case of chaotic synchronization related to the 
phenomenon of synchronization by noise [|r5,0 . 



iii) li p > P2 the Lyapunov exponent is positive, i.e. there exponential divergence of two 
initially close sequences C{t), behaving both chaotically. 

As a first check confirming that our advection-population dynamics model is close to the 
random map when mixing is strong, we fix tq = 1/2 and ri = 4, as in [|1^]. This describes a 
system in which predators are advected over regions in which not enough food is available 
{tq = 1/2 leads to population extinction) and over the strip-like regions in which preys are 
abundant (leading to chaotic population dynamics of the predators). We iterate (|5) to 
obtain C{i) = C(x(i),i), and then calculate the reaction Lyapunov exponent A^ for our 
system: 

1 N-l 1 N-l 

The initial condition C(x, 0) was a smooth function proportional to sin(x) sin(?/). A'^ mea- 
sures the rate of convergence or divergence of two initially similar concentration values C 
at the same fluid particle. In Fig. |l] we show A^ as a function of p for different values of K. 
When K is increased above a high enough value (e.g. K ^ 9, for which the standard map 
shows a unique ergodic chaotic region [jlll), A^ approaches the Lyapunov exponent A of the 
random map |12[. On the other side, when K is small, this correspondence is lost. 

Therefore, exploiting this equivalence we can study different regimes of our system de- 
pending on the value of the parameter p. In particular, the spatial patterns of the advected 
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field are strongly dependent of the value of A'^. Next section is dedicated to the study of 
these structures. 

IV. PREDATOR SPATIAL STRUCTURES 

The three regimes described above for the random map are also found for the behavior 
of the Lyapunov exponent A^ as a function of p in our advection model, with just some 
minor quantitative differences, e.g. in the values of pi and p2 (in particular, in most of our 
calculations we take K = 9, which gives pi ^ 0.34 and p2 — 0.48). These three regimes give 
rise to the following different predator spatial structures: 

i) For p < pi, the concentration of predators vanishes in all the space. The nutrient area 
is too small to support a stable population. 

ii) If pi < p < P2i a typical spatial pattern appears. A relatively high, but very intermit- 
tent, concentration of predators appears in the area occupied by preys. Moreover, the 
concentration pattern displays fractal features. 

iii) \i p > P2i the spatial concentration of predators shows a random pattern. No typical 
structure seems to emerge. 

Whereas the result for case i) is self-evident, cases ii) and iii) need a more elaborated 
study. We proceed in the following subsections. 

A. Case pi < p < P2 

The regime which we have labeled above with ii) is characterized by a negative reaction 
Lyapunov exponent A^. In this case, we observe numerically (see Fig. ^ the existence of 
a typical structure of the reactive field, which follows the strip-like structure of the preys. 
Moreover, a fractal pattern seems to be displayed by the distribution. Let us proceed to a 
quantitative characterization of these features. 
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References [jT5[ and |T6[ study the general continuous-time case of a chemically or bio- 



logically decaying field (thus with a negative reaction Lyapunov exponent) advected by a 
chaotic 2d flow. Since periodic velocity fields are used, it is straightforward to apply the re- 
sults in these papers to our case with discrete time. Nevertheless, a fundamental assumption 
in these studies is that a source term in the equations, analogous to our prey distribution, is 
a smooth function of space. The quantitative results of ||T5| , |l6| would fail (see e.g. Eq. (10) 



in ||T5|) when discontinuities are present in the source, as in our localized prey distribution 



(|13|). Therefore, in our calculations, and with the view on characterizing the spatial pat- 
terns, we will use a continuous approximation to the previous step function describing the 
distribution of nutrients, which would allow us to compare our results (in the regime of 



A < 0) with those in |T5|,|T6|. The approximation is performed by truncating the Fourier 



transform of the step function of width p (Eq. (|T3D and smoothing properly the coefficients 



1% . The final expression we use is: 



for any < x,y <27r and A^ = 10. Fig. (|^) shows a Id cut of this smoothed distribution of 
nutrients. 

A quantitative characterization of the observed structures can be performed in terms of 
structure functions. In particular, the structure function of order one, 5*1, is defined by: 

^i((5x) = (|C(x + 5x)-C(x)|), (18) 

where < ... > indicates an average taken over the different spatial points x along a line 



in the system. In [|T6| the scaling of Si is calculated with the result Si{6x.) ~ |(5x|" when 



5x -^ 0, with a ~ '-jjr when A^ < and A'^ > |A^|, being A^ the Lyapunov exponent of 
the Lagrangian motion (HI). The above expression for a is just an approximation to which 



multifractal corrections should be in principle added [|T^], but we are not going to consider 
them here. 

The Lyapunov exponent of the standard map, for K high enough, is given |]TI| by: 
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A^ ^ \n{K/2). In our calculations for i^ = 9 we are in the above mentioned conditions, 
that is, A-^ > |A-^| for all the values of p. 

Fig. d^) shows |A^| as a function of p for K = 9 (the smooth approximation to the 
spatial distribution of preys is used). In addition, we have numerically calculated the scaling 
exponent a of 5*1 in lines across the central strip of nutrients, and multiplied it by ln(9/2) ^ 
A^ for different values of p. The agreement between both quantities is quite good for p 
near pi, confirming the expression a ~ -^p^, although it gets worse as p ^ p2- The reason 
for this are the already mentioned multifractal corrections to the scaling exponent of the 
structure function, but a deeper discussion about this will be given in a subsequent work. 
The agreement allows us to understand the pattern displayed in Fig. (H) in terms of the 



filamental fractal patterns discussed in |15,16| for continuous-time dynamics. The observed 



fractal structures are revealing the stable and unstable manifolds (local contracting and 
expanding directions) attached to each point of the phase space of the standard map. It 
should be noted however that there is here a much larger amount of irregularities at small 
scales than in the patterns analyzed in [|I^,|r^. The reason is the much more irregular 



dynamics associated to the logistic map considered here. The Lagrangian evolution of C(t) = 
C(x(t), t) is related to the one of the logistic random map, which is of the on-off intermittency 
type. Models considered in |[15|JI6|1 displayed simple local relaxation behavior. The small- 



scale structure seen in Fig. ^ will introduce stronger multifractal corrections in higher order 
structure functions. We finally remark that when a discontinuous distribution of preys such 
as (|13D is considered, the relationship a ~ ^^ is not satisfied at all. 

B. Case p > P2 

Now we proceed to study the case p > P2- In this regime, the reaction Lyapunov exponent 
A^ is positive, i.e. the chemical or biological part of our system is also chaotic. 

The patterns calculated in this regime have random appearance, being dominated by 
strong small-scale irregularity with very small amount of structure (see Fig. |^). In fact, the 
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scaling exponents of the first-order structure function are close to zero, as corresponding to 
a random discontinuous field. 

This random structure is easy to understand once one has realized that A"^ > in this 
range of p: neighboring sites, even if they have initially nearly similar concentration values, 
and even when they remain close for long time so that they experience close values of the 
sequence x(t), will unavoidably develop growing differences in in concentration values, thus 
leading to the observed discontinuities at small scales. 

V. DISCUSSION 

Summing up, spatial structures with fractal features (of filamental type) appear for the 
predator field in a range of values of the size of the nutrient patch p. An increasing amount 
of small-scale randomness appears when p is increased, until structure is finally lost. The 
analogy with the random map model has allowed us to understand this behavior as being 
originated by the change in the value of the reaction Lyapunov exponent A^ when p is 
varied. In particular, structure is lost when A^ becomes positive. For p small enough, global 
extinction occurs, since most of the system has a parameter value for which C = is the 
only attractor. 

Our results have been obtained for a particular set of coupled maps, and for a specific 
nutrient distribution. We do not expect major qualitative changes in the above findings if 
the standard map is replaced by a different advecting flow, as long as the Lyapunov exponent 
A^ takes the same value. This belief is supported by the more detailed arguments of Refs. 



| T5| and [|1^] for the time-continuous case. It should be said however that the quantitative 



strength of multifractal corrections to simple expression such as a ~ -jjr will depend on the 
particular flow chosen. 

Our choice of the logistic map as the population dynamics to study is certainly important 
for the results obtained. Our results should describe the behavior under other population 
models as long as their parameters take values favoring chaotic oscillations in a localized 
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portion of space, and favoring relaxation to a fixed point in the rest. The election of the 
logistic map has allowed the use of results known for random logistic maps, thus helping 
to interpret the different patterns in terms of the value of A^ and its relationship with A^. 
Those quantities would be the right tool for the interpretation of advection-reaction patterns 
in other population or chemical models. 

Diffusion has been discarded in the present work. We expect that its only effect would 
be to smooth out any small-scale fractal or random structure below a size of the order of 



D/X^. In fact, our numerical calculations have an effective diffusion which comes from 
our minimal spatial resolution. As mentioned above, a more controlled way to introduce 
diffusion is to perform explicitly, after each map operation, an average of the concentrations 
of fluid particles closer than the diffusion length. 

We finally mention that the map approach turns out to be an extremely efficient method 
from the numerical point of view, as compared to direct solution of partial differential 
equations such as (|I]) or other continuous approaches. 
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Figure 1 
Lopez et al. 
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FIG. 1. Reaction Lyapunov exponent, A , calculated for different values of K. Solid line 

corresponds to K = 1.5, dotted line to i^ = 9 and dashed line to i^ = 19. 



Figure 2 
Lopez et al. 
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FIG. 2. Onedimensional cut of the continuous distribution of nutrients. 
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Figure 3 
Lopez et al. 
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FIG. 3. |A I for a continuous distribution of nutrients, calculated in two different ways against 
p. Values labeled with squares come from a direct calculation using expression (|16|) in the text. 
Circles are calculated from |A | ~ qA , being a the numerically calculated scaling exponent of the 
first order structure function. 
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Figure 4 
Lopez et al. 
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FIG. 4. Predator 2d pattern obtained for p = 0.37. The distribution of nutrients is continuous. 
The hghter the colour the higher the concentration. 
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Figure 5 
Lopez et al. 




1 2 3 4 5 6 

FIG. 5. Predator 2d pattern obtained for p = 0.87. The distribution of nutrients is continuous. 
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